---
title: "Expectation ratings analysis"
output: html_document
date: "2023-11-29"
---

```{r, include=FALSE}
library(base)       # The R Base Package
library(matrixStats)# Functions that Apply to Rows and Columns of Matrices (and to Vectors) 
library(sjmisc)     # Data and Variable Transformation Functions 
library(stringr)    # Simple, Consistent Wrappers for Common String Operations (filtering of partial name)
library(purrr)      # Functional Programming Tools

library(ggplot2)    # Create Elegant Data Visualisations Using the Grammar of Graphics
library(ggpubr)     # 'ggplot2' Based Publication Ready Plots
library(ggsignif)   # Significance Brackets for 'ggplot2'
library(ggsci)      # Scientific Journal and Sci-Fi Themed Color Palettes for 'ggplot2'
library(ggrain)     # A Rainclouds Geom for 'ggplot2'

library(psych)      # Procedures for Psychological, Psychometric, and Personality Research
library(irr)        # Various Coefficients of Interrater Reliability and Agreement
library(rstatix)    # Pipe-Friendly Framework for Basic Statistical Tests 
library(broom)      # Convert Statistical Objects into Tidy Tibbles 
library(lmerTest)   # Tests in Linear Mixed Effects Models
library(stats)      # The R Stats Package
library(lme4)       # Linear Mixed-Effects Models using 'Eigen' and S4 
library(emmeans)    # Estimated Marginal Means, aka Least-Squares Means
library(car)        # Companion to Applied Regression
library(pbkrtest)   # Parametric Bootstrap, Kenward-Roger and Satterthwaite Based Methods for Test in Mixed Models
library(effectsize) # Indices of Effect Size
```

```{r}
setwd("C:/Users/leu/OneDrive - UCL/Studies/Study - Expectation/Registered Report/Data acquisition/data/Expectation Matlab session file")
rat_exp_LL <- read.delim("rating_exp_low_low_all.csv", sep =",", header =FALSE)
rat_exp_HH <- read.delim("rating_exp_high_high_all.csv", sep =",", header =FALSE)
rat_exp_HM <- read.delim("rating_exp_high_med_all.csv", sep =",", header =FALSE)
rat_exp_LM <- read.delim("rating_exp_low_med_all.csv", sep =",", header =FALSE)

rat_perc_LL <- read.delim("rating_perc_low_low_all.csv", sep =",", header =FALSE)
rat_perc_HH <- read.delim("rating_perc_high_high_all.csv", sep =",", header =FALSE)
rat_perc_HM <- read.delim("rating_perc_high_med_all.csv", sep =",", header =FALSE)
rat_perc_LM <- read.delim("rating_perc_low_med_all.csv", sep =",", header =FALSE)

```

```{r}
rat_perc_LL <- stack(rat_perc_LL)
rat_perc_HH <- stack(rat_perc_HH)
rat_perc_HM <- stack(rat_perc_HM)
rat_perc_LM <- stack(rat_perc_LM)

means_HM <- describe(rat_perc_HM$values, ranges=TRUE)
means_LM <- describe(rat_perc_LM$values, ranges=TRUE)
rat_perc_all <- rbind(rat_perc_LL,rat_perc_HH,rat_perc_HM,rat_perc_LM)

#rat_perc_all$values_scaled <-10/1.524 *rat_perc_all$values #--> done directly in matlab to account for differences in VASmax
rat_perc_all$subject <- as.factor(rep(sprintf("Sub_%d",1:40),times=32))
rat_perc_all$trial <-as.factor(sort(as.numeric(rep(rep(sprintf("%d",1:8),times=40),decreasing=FALSE),times=4)))
rat_perc_all$cue <- as.factor(c(rep(("low"),times=320),rep(("high"),times=320),rep(("high"),times=320),rep(("low"),times=320)))
rat_perc_all$temperature <-as.factor(c(rep(("matched"),times=320),rep(("matched"),times=320),rep(("unmatched"),times=320),rep(("unmatched"),times=320)))


rat_exp_LL <- stack(rat_exp_LL)
rat_exp_HH <- stack(rat_exp_HH)
rat_exp_HM <- stack(rat_exp_HM)
rat_exp_LM <- stack(rat_exp_LM)
rat_exp_all <- rbind(rat_exp_LL,rat_exp_HH,rat_exp_HM,rat_exp_LM)

#rat_exp_all$values_scaled <-10/1.524 *rat_exp_all$values
rat_exp_all$subject <- as.factor(rep(sprintf("Sub_%d",1:40),times=32))
rat_exp_all$trial <-as.factor(sort(as.numeric(rep(rep(sprintf("%d",1:8),times=40),decreasing=FALSE),times=4)))
rat_exp_all$cue <- as.factor(c(rep(("low"),times=320),rep(("high"),times=320),rep(("high"),times=320),rep(("low"),times=320)))
rat_exp_all$temperature <- as.factor(c(rep(("matched"),times=320),rep(("matched"),times=320),rep(("unmatched"),times=320),rep(("unmatched"),times=320)))

```

```{r}
rat_exp_all$rating <- as.factor(rep("expectation"))
rat_perc_all$rating <- as.factor(rep("perception"))

merge_all <- rbind(rat_exp_all, rat_perc_all)
```

```{r}
wilcox.test(rat_exp_all$values)
wilcox.test(rat_perc_all$values)

check_data <- rat_perc_all
lmm <- lmer(values ~ cue*temperature+ (1|subject), data=check_data)
qqnorm(resid(lmm))
qqline(resid(lmm))

check_data$cook <- cooks.distance(lmm)
plot(check_data$cook, check_data$subject)
```


```{r}
ggboxplot (data=rat_exp_all, y="values", x="condition",
           size=1) 

ggboxplot (data=rat_perc_all, y="values", x="cue", group="cue", color="cue",facet.by = "temperature",
           size=1, add = "jitter") +
   geom_exec(geom_line, merge_all, x="cue", y="values", group="subject", color="darkgrey" )+
   labs(x="cue", y="ratings of intensity perception [VAS 0-10]")+
    font("legend.title", face = "bold")+
  font("xy.title", size = 16)+
  font("xy.text", size = 14)

 

ggboxplot (data=merge_all, y="values", x="rating",group="cue", color="cue", facet.by = "temperature",
           size=1, ylab="intensity rating VAS [0-10]")

ggplot(merge_all, aes(temperature, values, fill=cue, color=cue))+
  geom_rain(alpha=.5,rain.side = "r",
            boxplot.args = list(color = "black", outlier.shape = NA),
            boxplot.args.pos = list(position = ggpp::position_dodgenudge(x = .2), width = 0.15),
            violin.args.pos = list( side = "r",
              width = 0.7, position = position_nudge(x = 0.3)))+
    theme_classic() +
    scale_fill_brewer(palette = 'Set1')+
  labs(x="temperature", y="ratings of intensity perception [VAS 0-10]")+
  stat_compare_means(method = "t.test", paired=TRUE, label = "p.signif",label.x.npc = "middle")+
  facet_grid(cols=vars(rating))
#+
  stat_summary(fun = mean, geom = "line", aes(group = temperatur, color = temperatur)) +
  stat_summary(fun = mean, geom = "point", color="black",
               aes(group = rating))
  
  
```

ggplot(merge_all, aes(temperature, values_scaled, fill=cue, color=cue))+
  geom_rain(alpha=.5,rain.side = "r",
            boxplot.args = list(color = "black", outlier.shape = NA),
            boxplot.args.pos = list(position = ggpp::position_dodgenudge(x = .2), width = 0.15),
            violin.args.pos = list( side = "r",
            width = 0.7, position = position_nudge(x = 0.3)))+
    theme_classic() +
    theme(text = element_text(size = 20)) +
    scale_fill_brewer(palette = 'Set1')+
  labs(x="temperature", y="ratings of intensity perception [VAS 0-10]")+
  stat_compare_means(method = "t.test", paired=TRUE, label = "p.signif",label.x.npc = "middle")+
  facet_grid(cols=vars(rating))



ggplot(rat_perc_all, aes(temperature, values, fill=cue, color=cue))+
  geom_rain(alpha=.5,rain.side = "f2x2",
            boxplot.args = list(color = "black", outlier.shape = NA),
            boxplot.args.pos = list(position = ggpp::position_dodgenudge(x = 0.3), width = 0.15),
            violin.args.pos = list( side = "r",
              width = 0.7, position = position_nudge(x = 0.3)))+
    theme_classic() +
    scale_fill_brewer(palette = 'Set1')+
  labs(x="temperature", y="ratings of intensity perception [VAS 0-10]")+
  stat_compare_means(method = "anova", paired=TRUE, label = "p.signif")

plots <- rainplot_exp / rainplot_prec

ggplot(merge_all, aes(temperature, values, fill=cue, color=cue)+
  geom_rain(alpha=.5,rain.side = "r",
            boxplot.args = list(color = "black", outlier.shape = NA),
            boxplot.args.pos = list(position = ggpp::position_dodgenudge(x = .2), width = 0.15),
            violin.args.pos = list( side = "r",
              width = 0.7, position = position_nudge(x = 0.3)))+
    theme_classic() +
    scale_fill_brewer(palette = 'Set1')+
  labs(x="temperature", y="Visual Analog Scale [0-10]")+
  stat_compare_means(method = "t.test", paired=TRUE, label = "p.signif",label.x.npc = "middle")+
  facet_grid(cols=vars(rating))
  
### Linear Mixed Model for expectation ratings
```{r}
expect_lmm <- lmer(values~cue+(1|subject), data=rat_exp_all)
anova(expect_lmm, ddf="Kenward-Roger", type =3)
F_to_eta2( 6108 , 1, 1239) #cue
```

### Linear Mixed Model for perception ratings
```{r}
rating_lmm <- lmer(values~cue*temperature+(1|subject), data=rat_perc_all)
summary(rating_lmm)
if(requireNamespace("pbkrtest", quietly = TRUE))
#options(contrasts = c("contr.sum","contr.poly")) # sets contrast for R environment
anova(rating_lmm, ddf="Kenward-Roger", type =3)

emmeans(rating_lmm, list(pairwise ~ cue|temperature), adjust = "tukey" ) #post-hoc

F_to_eta2(3357.926, 1, 1237) #cue
F_to_eta2(  20.634, 1, 1237) #condition
F_to_eta2( 451.532, 1, 1237) #interaction


ggplot(merge_all, aes(temperature, values, fill=cue, color=cue))+
  geom_rain(alpha=.5,rain.side = "r",
            boxplot.args = list(color = "black", outlier.shape = NA),
            boxplot.args.pos = list(position = ggpp::position_dodgenudge(x = .2), width = 0.15), 
            violin.args.pos = list( side = "r",
              width = 0.7, position = position_nudge(x = 0.3)))+
    theme_classic() +
    scale_fill_brewer(palette = 'Set1')+
  labs(x="temperature", y="ratings of intensity perception [VAS 0-10]")+
  stat_compare_means(method = "t.test", paired=TRUE, label = "p.signif",label.x.npc = "middle")+
  facet_grid(cols=vars(rating))+
  font("legend.title", face = "bold")+
  font("xy.title", size = 16)+
  font("xy.text", size = 14)


```

ggplot(merge_all, aes(temperature, values, fill=cue, color=cue))+
  geom_rain(alpha=.5,rain.side = "r",
            boxplot.args = list(color = "black", outlier.shape = NA),
            boxplot.args.pos = list(position = ggpp::position_dodgenudge(x = .2), width = 0.15), 
            violin.args.pos = list( side = "r",
              width = 0.7, position = position_nudge(x = 0.3)))+
    theme_classic() +
    scale_fill_brewer(palette = 'Set1')+
  labs(x="temperature", y="Visual Analog Scale [0-10]")+
  stat_compare_means(method = "t.test", paired=TRUE, label = "p.signif",label.x.npc = "middle")+
  facet_grid(cols=vars(rating))+
  font("legend.title", face = "bold")+
  font("xy.title", size = 16)+
  font("xy.text", size = 14)

ggboxplot (data=rat_perc_all, y="values", x="cue", group="cue", color="cue",facet.by = "temperature",
           size=1, add = "jitter", legend="right") +
   geom_exec(geom_line, merge_all, x="cue", y="values", group="subject", color="darkgrey" )+
   labs(x="cue", y="intensity perception [VAS 0-10]")+
    font("legend.title", face = "bold")+
  font("xy.title", size = 16)+
  font("xy.text", size = 14)+
  theme(panel.border = element_blank(), axis.line = element_line()) +
    theme(axis.title.x=element_blank(), axis.text.x=element_blank(),axis.ticks.x=element_blank())
  